namespace hms::swe
{

template<typename MeshType>
template<BlockLocation loc>
void SolverThread<MeshType>::computeStorage(
	ArrayXXsMap& qBlock,
	ArrayXXsMap& hBlock,
	const Block<loc>& block,
	const Ref<const ArrayXXs>& dBlock,
	const Ref<const ArrayXXs>& zBlock,
	const Ref<const ArrayXXs>& vBlock
){
	for (Index y=0; y<block.nyRead(); ++y){
		using QMass = Eigen::Map<Array1Xs,0,Eigen::InnerStride<3>>;
		using QMom  = Eigen::Map<Array2Xs,0,Eigen::OuterStride<3>>;

		QMass   qMass { qBlock.col(y).data()  , 1, block.nxRead() };
		QMom    qMom  { qBlock.col(y).data()+1, 2, block.nxRead() };
		Array1XsMap h { hBlock.col(y).data()  , 1, block.nxRead() };

		Array1XsCMap d { block.template rowMapRead<1>(dBlock, y) };
		Array1XsCMap z { block.template rowMapRead<1>(zBlock, y) };
		Array2XsCMap v { block.template rowMapRead<2>(vBlock, y) };

		auto dRepl { d.template replicate<2,1>() };
		qMass = d;
		qMom  = (dRepl < SolverType::minWaterDepth).select( 0, v * dRepl );

		h = z + d;
	}
}

} // namespace hms::swe
